\subsection*{Jigsaw seeds and overlapping seeds}

With the heterogeneous seeds lookup table, the mapping costs of the seeds are
drastically reduced with little increase in memory overhead. However, this
solution is incomplete, as heterogeneous seeds require dividing the read into
fewer long seeds (at the length of $L_{max}$), which yields poor error tolerance
according to the pigeonhole principle.

If the read can be directly divided into cheap seeds at minimum lengths to query
the heterogeneous lookup table (rather than $L_{max}$), or if the long seeds can
overlap, as shown in Figure~\ref{fig:jigsaw} and in Figure~\ref{fig:overlap} respectively, then
we can still extract many seeds from the read, regaining the high error
tolerance. In this section, we propose ``jigsaw seeds" and ``overlapping seeds"
for the two methods respectively.

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{figures/JigSawSeeds_B.pdf}
\caption{An example of jigsaw seeds: a read is divided into cheap seeds with different lengths.}
\label{fig:jigsaw}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{figures/OverLappingSeeds_B.pdf}
\caption{An example of overlapping seeds: a read is divided into long cheap seeds overlapping each other.}
\label{fig:overlap}
\end{figure}

\subsubsection*{Jigsaw seeds: heterogeneous with higher error-tolerance}

The main difficulty of proposing a static seed partition policy with high
error-tolerance for heterogeneous seeds is that different seeds have different
lengths and their placements in the read is highly read dependent. If dividing
improperly, we may once again acquire expensive seeds. Luckily, with
heterogeneous lookup table, it is unnecessary to have a static seed partition
policy, instead, we propose a cheap dynamic seed partition policy that
explores the minimal lengths of the cheap seeds to obtain the maximum number of
seeds, which offers higher error-tolerance.

Instead of dividing the read into seeds in parallel, we propose a dynamic policy
which divides the read into seeds in serial, as shown in Algorithm~\ref{alg:SeedPar}. While
each time the mapper still uses a long seed ($L_{max}$ bases) to query the
heterogeneous lookup table, with the exact matching extension length returned by
the lookup table, the mapper can reduce the length of the querying long seed
($L_{max}$ bases) to the minimal length for the seed to remain cheap. The rest
of the bases of the long seed is returned to the read to form the next seed. In
this way, a read is partioned into as many non-overlapping cheap seeds as
possible.

%\begin{algorithm}
%\caption{A dynamic seed partitioning algorithm which partitions the read into as
%many non-overlapping seeds as possible}
%\label{alg:SeedPar}
%\begin{algorithmic}[1]
%\Procedure{JigsawSeeds}{$READ$}\Comment{Partition READ into SEEDs}
%	\State $i \gets 1$
%	\State $START \gets 0$\Comment{Set the head of the first SEED to 0}
%	\While{$i \leq e + 1$}
%		\State $SEED \gets READ.substr(0, L_{max})$\Comment{Get the
%			first SEED with length $L_{max}$}
%		\State $Length \gets HET_LOOKUP_TABLE.query(SEED)$\Comment{Get
%			the exact length of the matching extension of the seed}
%		\State $SEED \gets SEED.substr(0, LENGTH)$\Comment{Recalibrate
%			the length of SEED to LENGTH}
%		\State $START \gets START + LENGTH$\Comment{Set the head of the
%			next SEED}
%	\EndWhile
%\EndProcedure
%\end{algorithmic}
%\end{algorithm}

%\mbox{}\\

\begin{algorithm}[h]
\DontPrintSemicolon
\KwData{HET\_LOOKUP\_TABLE}
\KwIn{READ} 
\KwResult{Obtain $e + 1$ SEED}
 $i=0$\
 $START=0$\
 \While{$i < e + 1$} {
 	SEED = READ.substr(START, $L_{max}$)	\tcp*[l]{Obtain the long seed with the length $L_{max}$ at the START position}
\mbox{}\\
	LENGTH = HET\_LOOKUP\_TABLE.query(SEED)	\tcp*[l]{Get the exact length of the matching extension of the long seed}
\mbox{}\\
	SEED = SEED.substr(0, LENGTH)	\tcp*[l]{Recalibrate the length of the seed to match the extended cheap long seed}
\mbox{}\\
	START += LENGTH	\tcp*[l]{Increment the START position of the next long seed by LENGTH bases}

 }
\caption{A dynamic seed partitioning algorithm that partitions the read into as
many non-overlapping seeds as possible}
\label{alg:SeedPar}
\end{algorithm}

%\mbox{}\\

With the above algorithm, the mapper can obtain the maximum number of cheap
seeds with different lengths jostling each other in the read. This is akin to a
jigsaw puzzle, where each piece of the puzzle is all different but they
interlock with each other, building up the entire picture. Therefore, we call
this dynamic seed partition policy ``Jigsaw Seeds".

Unlike hobbes~\cite{hobbes}, which needs many accesses to the lookup table in order
to find the set of cheapest seeds, jigsaw seeds only need $N$ accesses for $N$ seeds
and they guarantee that all seeds are as cheap as possible (before hitting the
extension limit).

%For each read, we use the first few bases to form a virtual
%long seed to query the lookup table, which returns the effective length of the
%seed together with the locations of the seed. Then we crop the virtual long seed
%to its effective length, which yields the jigsaw seed.  Right at the end of the
%first jigsaw seed, we continue with  the second virtual long seed and use it to
%query the lookup table as well, until we cover the entire read with $e+1$ jigsaw
%seeds to tolerate $e$ errors.

\subsubsection*{Overlapping seeds: a backup when jigsaw seeds fails}

Jigsaw seeds works most of the time, since usually there are few long cheap
seeds in a read. Nonetheless, there are still a few cases where a read contains
multiple long seeds. For those reads, even partitioned with jigsaw seeds,
there are still not enough seeds to tolerate all possible errors.

The problem originates from the fact that the pigeonhole principle requires the
seeds to be non-overlapping. If the seeds are allowed to overlap each other
slightly, where at any given point in the read, there are at most two seeds
overlapping each other (Figure~\ref{fig:overlap}), then we can still partition
the read into many seeds, even if there are multiple long cheap seeds. We call
such partition policy as \textbf{overlapping seeds}. However, it is now
possible that a single error simultaneously destroys seeds at once. As a
result, $e$ errors may destroy all $e+1$ overlapping seeds. Fortunately, as we
will describe below, we can provide a generalization of the pigeonhole principle
by allowing seeds to have a single-base error in the extension part,
restoring the high error-tolerance property of the mapper.

\begin{overlap} \label{overlap}

If a read contains $e$ errors and is partitioned into $e+1$ overlapping seeds,
where there are at most two seeds overlapping each other at any point in the
read, then there must be at least one seed with no errors in its non-overlapping
region and at most one error in its overlapping region.

\end{overlap}

\begin{proof} 

% CALKAN: stopped here on Sat, Aug 17, 19:03 EEST.
We prove by contradiction. Let us assume that such seed does not exist. Then
each seed will have at least one error in its non-overlapping region or have
at least two errors in its overlapping region. For the first case, if one or
more errors is in the non-overlapping region of a seed, then no other seed
shares the error, putting at least one error per seed. For the second
case, since there are exactly two seeds sharing a single overlapping region,
each seed inherits one half of the error. Consequently, at least two errors in
the overlapping regions of two seeds is also equivalent to at least one error
per seed of the two seeds. As a result, $e+1$ seeds will have at least $e+1$
errors, but we have only $e$ errors, which is a contradiction.

\end{proof}

\begin{overlapLemma} 

With the same conditions as above, there must be a seed which not only has
at most a single error in the overlapping region, but also the error
locates specifically in its right overlapping region. In other words, there
exists a seed which has no error in its left-overlapping region (if the seed
overlaps with its left neighbor seed) or in its non-overlapping region.
\label{seedlemma}
\end{overlapLemma}

The proof of the Lemma is provided in the Appendix.

We further observe that in the heterogeneous lookup table, most expensive seeds
have only a few extensions (12 on average with 3124 as the maximum). This
observation suggests that instead of searching for perfect extensions of the
long seed in the extension table, it is computationally inexpensive to search
for extensions while allowing one error in the extension part of the seed.
Although allowing errors forces the mapper to give up cheap searching algorithms
(i.e., binary search), the additional computation is still much cheaper than
simply using the short seeds. For example, for the expensive short seed that has
the most extensions (3124), a linear search with one allowed error generates
3124 extension verifications (a small-scale verification for only the extension
part of the seed), which returns one (or sometimes a few) cheap long seed with
fewer than 100 locations. Compared using short seeds, where the same expensive
short seed directly returns 244891486 locations to verify, the total
computational cost is much less.

Moreover, unlike the verification calculation for the entire read, the
extension verification need not to be 100\% correct, since any extension
passes the small-scale verification will go through a full read-long-scale
rigorous verification again later on. Extension verification only serves the
purpose of a filtering mechanism. As long as the small-scale
verification does not falsely reject correct extensions (as false negatives), it
is acceptable for the verification to allow a few incorrect extensions (with
more than one error) pass (as false positives).

We also propose a simple vectorization-friendly algorithm for extension verification,
as shown in Algorithm~\ref{alg:extVerify},
which runs very fast on modern processors that are capable of vector processing.

%\clearpage
%\begin{algorithm}
%\caption{A vectorization-friendly extension verification algorithm which verifies if
%	there exists more than 1 error in the extension part of the seed}
%\label{alg:extVerify}
%\begin{algorithmic}[1]
%\Procedure{ExtensionVerify}{$SeedExt TableExts
%	ACTATCGT$}\Comment{Verify if SeedExt contains at most 1 error from
%		TableExt}
%	\State$hamming\_mask \gets SeedExt \oplus TableExt$\Comment{Calculate
%		the bit-mask of matching bases;}
%	
%	\State$insertion\_mask \gets (SeedExt \oplus (TableExt \gg 1)) \land hamming\_mask$\Comment{Calculate
%		the bit-mask of matching bases with 1 insertion;}
%		
%	\State$hamming\_mask \gets SeedExt \oplus TableExt$\Comment{Calculate
%		the bit-mask of matching bases;}
%	
%		\State $START \gets 0$\Comment{Set the head of the first SEED to 0}
%	\While{$i \leq e + 1$}
%		\State $SEED \gets READ.substr(0, L_{max})$\Comment{Get the
%			first SEED with length $L_{max}$}
%		\State $Length \gets HET_LOOKUP_TABLE.query(SEED)$\Comment{Get
%			the exact length of the matching extension of the seed}
%		\State $SEED \gets SEED.substr(0, LENGTH)$\Comment{Recalibrate
%			the length of SEED to LENGTH}
%		\State $START \gets START + LENGTH$\Comment{Set the head of the
%			next SEED}
%	\EndWhile
%\EndProcedure
%\end{algorithmic}
%\end{algorithm}


%\begin{center}
%\line(1,0){460}
%\end{center}

\mbox{}

\begin{algorithm}[H]
\DontPrintSemicolon
\KwIn{Two short sequences $SeedExt, TableExt$	\tcp*[h]{As a running example, assume
       SeedExt=ACTACGTT, TableExt=ACTATCGT} }
\KwOut{A boolean if the edit
distance is smaller or equal to 1 between the two sequences}

\mbox{}\\

$hamming\_mask \gets SeedExt \oplus TableExt$\tcp*[h]{Compute the bit-mask
for the matching bases. For the running example, hamming\_mask = 00001111; 0
represents a matching base while 1 represents a mismatching base}

\mbox{}\\

$insertion\_mask \gets (SeedExt \oplus (TableExt \gg 1)) \land hamming\_mask$
\tcp*[h]{ Compute the bit-mask for the matching bases with 1 insertion. For the
running example insertion\_mask = 00001111; 0 represents a matching base while
1 represents a mismatching base or an inserted base}

\mbox{}\\

$deletion\_mask \gets ((SeedExt \gg 1) \oplus TableExt) \land hamming\_mask$
\tcp*[h]{Compute the bit-mask for the matching bases with 1 deletion. For the
running example
deletion\_mask = 00001000; 0 represents a matching base while 1 represents a
mismatching base or a deleted base}

\mbox{}\\

\tcp{Counting errors}

\mbox{}\\

\uIf{$count\_ones(hamming\_mask) \leq 1$}{

	\Return{TRUE}	\tcp*[h]{Return true if there is none or only one
mismatch.}
}

\uElseIf{$count\_ones(insertion\_mask) \leq 1$}{
	\Return{TRUE}	\tcp*[h]{Return true if there is none or only one
deletion.} 
}

\uElseIf{$count\_ones(deletion\_mask) \leq 1$}{
	\Return{TRUE}	\tcp*[h]{Return true if there is none or only one
insertion. For the running example, TRUE is returned.} 
}

\uElse{
	\Return{FALSE}	\tcp*[h]{Return false as there are more than one errors} 
}

\caption{A vectorization-friendly extension verification algorithm that verifies if
	there exists at most 1 error in the extension part of the seed}
\label{alg:extVerify}
\end{algorithm}

\mbox{}\\

Notice that the above algorithm allows false negatives. For example, AAAACCAAAA
is regarded containing one error from AAAAATAAAA while in fact it contains two.

With the above observations, we complement the overlapping seed partition policy with
a requirement of a minimum 12 bases (or the same length of the short
seeds in the heterogeneous lookup table) non-overlapping and left-overlapping joint
region for each overlapping seed. This additional requirement guarantees the
existence of a seed which contains no error in its basic part and at most one
error in its extension part.

With overlapping seeds, we approximate the error-threshold of short seeds with
heterogeneous seeds.

 

%
%for each seed
%overlapping with other seeds, we require the total length of the left
%overlapping region and the non-overlapping region of the seed exceeds the length
%of the short seed length (12 bases) of the heterogeneous lookup table. We call
%this 
%
%
%
%
%****************Hongyi: Start revising portion:
%
%With the generalized Pigeonhole Principle and the Lemma above, we can redivide
%the read as the following: we first divide the read into ``basic seeds" that are
%as long as the basic seed length in the heterogeneous lookup table. Then we
%extend each ``basic seed" to its right to make them overlapping long seeds.
%Notice that the extension part of each long seed overlaps with the basic part of
%its right neighbor seed as shown in Figure~\ref{fig:overlap}. Recall earlier when we
%generate the heterogeneous seed lookup table, we set the maximum extension
%length as one base from the basic seed length. That is to guarantee that no seed
%can extend beyond the basic part of its right neighbor, which enforces the
%requirement that there can be at most two seeds overlapping each other at a
%time.
%
%If there are $e$ variations and we have $e+1$ overlapping long seeds, as the
%lemma suggested, we know that there must exist a seed which has no error in the
%basic part (the combination of left overlapping region and the non-overlapping
%region) and has at most one error in the extension part (the right overlapping
%region), no matter how the errors are distributed. As a result, for each
%overlapping seed, if we can find all locations where the reference sequence at
%the location has at most one error from the long seed in the extension part, we
%guarantee that we will not miss any of the correct mapping of the read. This is
%equivalent to finding all seeds in the lookup table, which contains at most a
%single error from the target seed in the extension part.  Luckily, the
%heterogeneous lookup table we proposed earlier fulfills this requirement.
%
%When we query the heterogeneous lookup table, since it contains no error in the
%basic part of the seed, we can use the it as an index for a direct access to
%find the basic seed in the lookup table.  For the extension part, however, as we
%need to tolerate one potential error, we cannot just find a perfect match, but
%rather we have to find all potential extensions in the extension table that are
%one edit away from the extension part of the target seed. The simplest way is to
%do a linear search within the extension table which compares every key, which is
%one of the possible extensions of the basic seed in the reference genome,
%against the extension part of the target seed. Fortunately, most of the basic
%seeds do not have a lot of extensions and for the few seeds that do have many
%extensions (usually are fairly frequent in the genome as well), we also propose
%a fast searching algorithm, which utilizes the space-time tradeoff, in the
%Attachment.
%
%Even with linear search, overlapping seeds saves a lot of computation.  Each
%extension comparison, which is merely a very small-scale edit distance
%calculation that only compares around ten more bases, eliminates tens to
%hundreds of locations which avoids tens to hundreds of full read-long
%edit-distance calculations.
%
%Furthermore, when there can be at most only one error, instead of using
%Smith-Waterman edit-distance algorithm, we can simply use a base-wise comparison
%with shifting algorithm to check if there are more than one edit. The
%pesudo-code is provided below:
%
%\begin{algorithm} \DontPrintSemicolon \KwIn{Two short sequences $\langle SeedExt
%\gets ACTACGTT, TableExt \gets ACTATCGT \rangle$} \KwOut{A boolean if the edit
%distance is smaller or equal to 1 between the two sequences} \tcp{Finding
%matches between sequences, 0 is a match, 1 is a mismatch} $hamming\_mask \gets
%SeedExt \oplus TableExt$\; \tcp{hamming\_mask = 00001111} $insertion\_mask \gets
%(SeedExt \oplus (TableExt \gg 1)) \land hamming\_mask$\; \tcp{insertion\_mask =
%00001111} $deletion\_mask \gets ((SeedExt \gg 1) \oplus TableExt) \land
%hamming\_mask$\; \tcp{deletion\_mask = 00001000} \tcp{Counting mismatches}
%\uIf{$count\_ones(hamming\_mask) \leq 1$}{ \tcp{The two sequences have one or
%zero substitution} \Return{TRUE}\; } \uElseIf{$count\_ones(insertion\_mask) \leq
%1$}{ \tcp{The two sequences have one or zero insertion} \Return{TRUE}\; }
%\uElseIf{$count\_ones(deletion\_mask) \leq 1$}{ \tcp{The two sequences have one
%or zero deletion} \Return{TRUE}\; } \uElse{ \tcp{There are more than one edits}
%\Return{FALSE}\; } \end{algorithm}
%
%With the help of modern processor which equipped with vector processing power,
%the above base-wise comparison can be done in just a few processor cycles.
%
%The maximum number of allowed errors $e$ of a read depends on how many
%``basic seeds" the read can be divided to. If a read can be divided into $N$
%``basic seeds", then as we proved above, the read can contain at most $N-1$ errors
%and will still be mapped correctly. In the heterogeneous lookup table, ``basic
%seeds" are essentially short seeds, thus the error tolerance of
%overlapping long seeds at least equals the error tolerance of short seeds.
%
%With heterogeneous lookup table, jigsaw seeds and overlapping seeds, we
%successfully attempted the benefit of long seeds, which provide fewer potential
%locations to verify, while retained the advantage of short seeds---without compromising the error tolerance of the mapper
%or incurring a lot of memory overhead.
%
%
%******************Hongyi: End revising portion.
